Nonlinear structures and thermodynamic instabilities in a one-dimensional lattice 

system 
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The equilibrium states of the discrete Peyrard-Bishop Hamiltonian with one end fixed are com- 
puted exactly from the two-dimensional nonlinear Morse map. These exact nonlinear structures 
are interpreted as domain walls (DW), interpolating between bound and unbound segments of the 
chain. The free energy of the DWs is calculated to leading order beyond the Gaussian approxi- 
mation. Thermodynamic instabilities (e.g. DNA unzipping and/or thermal denaturation) can be 
understood in terms of DW formation. 

PACS numbers: 87.10.-|-e, 63.70.-|-h, 05.70.Jk, 05.45-a 



Computational results 0] have suggested, at varying 
levels of numerical rigor, that a class of one-dimensional 
lattice models proposed within a variety of physical con- 
texts (e.g. interfacial wetting Q or DNA denaturation 
0) may exhibit thermodynamic instabilities Q which 
share many of the properties of ordinary phase transi- 
tions. 

Owing to the lack of analytical results, it is important 
to explore alternative strategies which reveal exact fea- 
tures of these instabilities. One such strategy, which will 
be pursued in this work, is to exploit some of the tools 
of nonlinear dynamics in order to study the properties of 
the underlying exact equilibrium structures. 

The general class of models under consideration is de- 
scribed by a Hamiltonian with a configurational part 
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where ?/„ is the transverse displacement of the nth site, 
i? is a dimensionless coupling constant and V{y) is any 
potential with a repulsive core, a stable minimum and a 
flat top. We will deal with the case of a Morse potential 
0, V{y) = (1 — e^'^Y . All quantities referred to in this 
paper are dimensionless. 

We will demonstrate the existence of exact nonlin- 
ear structures which correspond to domain walls (DW)s, 
"interpolations" from the bound to the unbound phase. 
It will further be shown that the detailed stability proper- 
ties of the DWs essentially determine the system's behav- 
ior with respect to both mechanical and thermal instabil- 
ities (unzipping and thermal melting, respectively, in the 
DNA context). Finally, we will use thermodynamic per- 
turbation theory to account for the minor discrepancies 
which occur between the predictions of Gaussian DW- 
based theory and standard numerical results based on 
the transfer integral (TI) method. This work describes 
the case i? ^ 1 (extreme discretization) ; DWs in the case 
i? <C 1 (continuum limit) have been treated in Ref. I^. 
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FIG. 1: The unstable manifold of the FP of the map ^ 
for R = 10.1 and A'' = 28. Black squares belong to stable 
equilibria, red open circles belong to unstable equilibria. The 
horizontal line at j/ = 44 demostrates the multivaluedness 
of the manifold as a function of y (4 stable and 3 unstable 
equilibria with that value of y; details in the inset). The 
vertical lines are drawn at pmin and Pmax, the minimal and 
maximal asymptotic slopes of DWs. 

The equilibria {yi"''} of J^l, subject to fixed-end 
boundary conditions yoj2/JV+i — const, are the solutions 
of the two-point boundary value problem for the second 
order recurrence relation yn+i — 2yn + y-n-i = RV'{yn)- 
By introducing p„ = yn—yn-i-, this is equivalent to orbits 
of the two-dimensional map 
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where n — 1, ■ ■ ■ , N , yi — pi + yo, and pi is unspecified, 
i.e. it can take any value which leads to the given yN+i 
after N iterations. The map Q has a single hyperbolic 
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FIG. 2: The 8 stable equilibria corresponding to = 28, 
1/0 = 0, yN+i = L = 80. Not shown are 7 unstable equilibria 
enmeshed between the stable ones. Inset: total energies for 
both stable (black squares) and unstable (red open circles) 
equilibria. The continuous curve corresponds to a theoreti- 
cal estimate which does not distinguish between stable and 
unstable equilibria (cf. text). Note that the typical energy 
difference - which does not shrink with increasing L - between 
neighboring extrema is about 0.15. 



fixed point (FP) at {p 
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0). Equilibria can be 



classified as stable or unstable, according to whether all 
eigenvalues of the Hessian 



dyndy„ 



(3) 



are nonnegative or not. 

The tangent map at the FP is given by 
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its eigenvalues are A± = 1 + R zt ^/B? + 2R and 
its corresponding (unnormahzed) eigenvectors k^^^ = 
(2_R, R ± \/-R2 + 2i?) , where the upper (lower) sign corre- 
sponds to the unstable (stable) manifold. The structure 
of the unstable manifold of the FP is shown in Fig. ^ 
Its multivaluedness as a function of y turns out to have 
important consequences when the fixed-ends boundary 
conditions appropriate to the problem are introduced. 
These are of two types: (a) yo — yN+i — 0, and (b) 
yo = 0,yN+i = L- The only equilibrium compatible 
with (a) is the FP. In order to construct the latter, we 
generated a large number of sequences using the initial 
conditions yi = pi = kJ/k^ exp(— s) with s uniformly 
spaced in an interval (si, S2) and initially accepted a se- 
quence if I^AT+i — L\/L < 2 • 10~^; a further scan of ac- 
cepted sequences was made to eliminate sequences which 
are "immediate neighbors" in the manifold. The result 
for L = 80 is a total of 8 stable and 7 unstable such 
sequences. They have the form of DWs consisting of a 
segment bound in the well of the Morse potential con- 
nected to a free segment on its plateau. The stable ones 



are shown in Fig. |2| at the left they more or less coincide 
with the fixed point = 0; to the right, they emerge as 
a bundle with linearly growing displacements and a range 
of slopes < Pa < Pmax- The Corresponding energies, 
according to 1^ are shown in the inset as functions of the 
final slope p. Note that the unstable equilibria have ener- 
gies slightly above the neighboring stable ones. The over- 
all dependence of the energy on the slope p can be under- 
stood by the following simple argument: each unbound 
site contributes an amount l+p"^ / {2R) to the energy; for 
L 3> 1, the number of unbound sites is approximately 
equal to L/p; hence E{p,L) = {l/p + p/2R)L; this is the 
dotted curve shown in the inset of Fig. |21 It follows that 
the minimum energy occurs at p — p* = {2R)^/^ and has 
value E*{L) = 2L/p*. 

A few remarks are in order here: (i) The number of 
equilibria grows linearly with L; more precisely, for suffi- 
ciently large L, the number of local minima is [L/pmin] — 
[L/pmax] + 1; there is one fewer unstable equilibrium, (ii) 
There are no zero eigenvalues for any equilibria, except 
for the discrete values of L at which a new pair is cre- 
ated, (iii) The unstable equilibria have only one nega- 
tive eigenvalue, so we call them saddles; the rest of their 
spectrum is essentially identical with that of neighboring 
local minima, (iv) For large R, the limits Pmin,Pmax are 
asymptotically Pmin ~ ln(2i?) -I- 1, pmax ~ + In 2, 
which we derived by approximate solution of the equa- 
tions for the unstable manifold, (v) The energy spacing 
(Peierls barrier) between the absolute minimum and the 
saddle with the lowest energy (cf. inset in Fig. |2Jl is 
about 0.15. (vi) For any given equilibrium (DW) with 
a transverse displacement L 3> 1 and final slope p, the 
quantity dE/dL = V {L) + {L — yt^)/ R p/R represents 
the force acting on the end particle. In the case of the ab- 
solute minimum, p — p* , this is equal to /o = (2/i?)^/^; 
in the context of DNA denaturation models, it corre- 
sponds to the force required to bring about mechanical 
"unzipping" at zero temperature. 

What happens at nonzero temperatures? In order to 
examine the relevance of DWs for thermodynamics, it is 
necessary to look at the partition functions .^Ar(O) and 
Zjv(L) which correspond to the two types of boundary 
conditions above. The quantity of interest is the differ- 
ence in the free energies 
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in particular, the derivative {dAG/dL)rp is equal to the 
force f{L,T) required to maintain the right end of the 
chain at a given tranverse displacement L. For small 
displacements around any equilibrium, 2/„ = 
$ can be expanded as 
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FIG. 3: Eigenvalue spectra of the hessians A^'^'' (open 
squares) and A^*^ (open circles) for A'^ = 32 and L = 100. 
The DW's spectrum consists of bands of optical and acous- 
tic phonons, localized respectively in the bound and unbound 
portions of the chain, and a single local mode in the gap; 
both bands are well described (to order 0(1/1/)) by the cor- 
responding free phonon dispersion curves (dotted). 



where Vm]a = V^^\ym'')/j\, and terms of higher than 
fourth order have been dropped. The partition function 
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can be calculated in the Gaussian approximation by keep- 
ing terms up to second order in l^. The result is 
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where {A^f^}, v = 1,2,..., TV are the eigenvalues of the 
Hessian matrix A'"^ . 

Eq. ^ will be apphed twice. First, to evaluate Zjv(O); 
in this case, the relevant equilibrium is the FP. Second, 
to evaluate Zn{L), in which case the relevant equilibrium 
is the DW with the minimal energy and the final slope 
p* (cf. above). Note that at temperatures of order 1 or 
higher (relevant for the occurrence of macroscopic insta- 
bilities), the roughness of the energy landscape (Fig. [SJ, 
typified by the height of the Peierls barrier (cf. comment 
(v) above), essentially disappears. 

The ratio of partition functions in jSJ can be approxi- 
mated as 
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The eigenvalue spectra { A^,*-* , aI,"^ } are shown in Fig. El 
As has already been pointed out in the context of the 
continuum-limit|^, the DW acts as an interface between 
bound and imbound phases; the unbound phase supports 
only acoustic phonons and the bound phase only optical 



ones. The quantity 

f) = i^ln{A(°)/AW} = iicT + 0(l), (10) 



where cr = In 
L/p* > 1 and iV 



(^VW2+ + ^/2)> is found to be, for 
L/p* ^ 1, independent of N and 
proportional to the number L/p* of unbound sites; the 
reason is that, as can be seen from Fig. |3| the optical 
phonons in the bound region of the DW cover exactly 
the frequency range of the optical phonons around the 
fixed point, and therefore cancel out in taking the ra- 
tio; this leaves only a dependence on the ratio between 
acoustic and optical frequencies for a number of modes 
of order L/p*. The proportionality constant cr(i?) has 
been calculated in Ref. p using the argument sketched 
above, without reference to the exact eigenvalue spec- 
trum. Irrelevant corrections of order unity originate from 
the transition region of the DW and/or the single local- 
ized oscillation mode. Inserting H1U|I and in Q, we 
obtain, in the Gaussian approximation. 
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It follows that the force now required to maintain the 
chain at (any) transverse displacement L 3> p* , 
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p* 
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is reduced as a result of thermal motion. At T = 2/cr this 
generalized unzipping force vanishes; spontaneous forma- 
tion of the DW becomes possible. In the context of DNA 
denaturation theory, this amounts to spontaneous melt- 
ing (thermal denaturation). That transition has been 
extensively studied, assuming periodic boundary condi- 
tions, both in the continuum and the discrete regime by 
means of the TI method. The critical temperature 2/ a 
obtained in the Gaussian approximation is in good agree- 
ment with the TI results for values of ii < 4; systematic 
deviations occur at higher R however 

It is possible to improve the calculation of AG by going 
beyond the Gaussian approximation. Standard thermo- 
dynamic perturbation theory 6] can be applied to eval- 
uate the contribution of cubic and quartic terms in © 
to the partition function ||7J), - again, both for the FP 
and the DW - and to derive a low-temperature expan- 
sion for the resulting free energy. Each power of the 
displacement field in (jSJ contributes a factor T^/^. The 

irreducible graphs which contribute to In Z^^^ to leading 
order in the temperature are: 

0«0 = -3r^«Wj5^^}2 ^and (13) 
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FIG. 4: Lowest order anharmonic correction to the free 
energy difference AG: The coefficient Ac2(p, L) for L = 
100, 140, 200, 300, 450 vs. p; for each L the difference (C2(p) - 
C2)/L is numerically evaluated at all stable minima; the in- 
terpolates taken at p = p* are plotted against 1/L (inset); 
this procedure extrapolates to a value Ac^ = —0.01784(5). 




unbound phase 



FIG. 5: Phase diagram in the ,f — T plane. The dotted 
line is obtained on the basis of the Gaussian approximation 
(predicted Tc = 1.291). The solid line is Eq. [TSlwith Acj = 
-0.01784 (cf. Fig. BJ. Inset: detail of critical region, showing 
a predicted Tc = 1.215; the triangle is the result of a finite-size 
scaling numerical TI calculation Tc — 1.227 



where gmX is the relevant Green's function (inverse of the 
Hessian evaluated at the equihbrium a). 

Let —TC2 denote the sum of the two contributions 
above. The correction to the free energy AG will be 
obtained as a difference between values at the two rele- 
vant equilibria (FP and DW with minimal energy); the 

quantity — is expected to be independent of 
N and proportional to the length of the unbound seg- 
ment (cf. above). Fig. ^ (inset) displays the values 

of Ac2(35*) = {C^*^ - 0^2^) /L (interpolated from the 
values at the closest DW's, displayed in the main fig- 
ure) for a series of L-values. The extrapolated value 
Ac2 = lim^^oo Ac2(p*, -/j) can be used as a correction 
to eqs. II II and [T^ above, i.e. 
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The resulting phase diagram is shown in Fig. |S1 The 
predicted — 1.215 is in good agreement with the nu- 
merical TI calculation [J;] TJ' = 1.227. 

We conclude with two remarks. The first is a caution 
about the inherent limitations of the low-order perturba- 
tional approach used. At lower values of i?, which tend 
to increase Tc to values substantially higher than 1, it 
cannot be expected to work (although the Gaussian ap- 
proximation produces very good estimates of Tc). At 
much higher values of i?, the perturbation terms in the 
Hamiltonian become too large compared with the Gaus- 
sian approximation; again, one would presumably have to 
include higher order terms to obtain reasonable results. 

The second remark concerns a further prospect: in- 
cluding nonlinear stacking terms in the Hamiltonian 
has been known to produce real or apparent first order 
transitions. It would be interesting to explore how this 
behavior would be reflected in the DW approach. 
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contract HPRN-CT-1999-00163 (LOCNET network) is 
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